/*
 * testComponents.cpp
 *
 *  Created on: Mar 30, 2011
 *      Author: curt
 */

#include "CUDA_vec.h"
#include "CUDA_full_weighting.h"
#include "CUDA_interpolation.h"
#include "CUDA_JacobiRelax.h"
#include "CUDA_Poisson.h"
#include "ConjugateGradient.h"
#include "utility.h"
#include <math.h>
#include <stdio.h>
typedef CUDA_Full_Weighting RestrictionOp;
typedef CUDA_Interpolation InterpolationOp;
typedef CUDA_JacobiRelax RelaxOp;
typedef CUDA_Poisson Matrix;
typedef CUDA_vec<float> Vector;

float epsilon = 1e-3;

int size(int p) { return (1<<p) - 1; }

bool test_norm(int p)
{
	bool passed = false;
	int N = size(p);
	Vector b(N*N);
	b.set_constant(1);
	float nm = norm(b);
	passed = abs(nm - N) < epsilon;
	return passed;
}

bool test_restrict(int p){
	bool passed = false;
	int N = size(p);
	int Nc = size(p-1);
	float * test_array = new float[N*N];
	for(int i=0; i<N; i++)
		for(int j=0; j<N; j++)
			test_array[i * N + j] = ((i+1)&0x1 != 0) || ((j+1)&0x1 !=0);

	Vector fine_x(N*N);
	Vector coarse_x(Nc*Nc);

	fine_x.copyToDevice(test_array, N*N);
	RestrictionOp restrict_op;

	restrict_op(fine_x, coarse_x);
	float nm = norm(coarse_x);
	passed = abs(nm - 0.75 * Nc) < epsilon;
	printf("Restricted norm: %f Expected: %f\n", nm, 0.75 * Nc);
	return passed;
};

bool test_interpolate(int p){
	bool passed = false;
	int Nc = size(p);
	int N = size(p+1);
	Vector coarse_x(Nc*Nc);
	Vector fine_x(N*N);

	coarse_x.set_constant(1);
	InterpolationOp interpolate;

	interpolate(coarse_x, fine_x);
	float nm = norm(fine_x);
	passed = abs(nm - N) < epsilon;
	printf("Interpolated norm: %f Expected: %f\n", nm, (float)N);
	return passed;
}

bool test_matrix(int p){
	bool passed = false;
	int N = size(p);
	Vector x(N*N);
	Vector Ax(N*N);
	x.set_constant(1);
	Matrix A;
	A(x,Ax);

	float nm = norm(Ax);
	passed = abs(nm - 16.1245) < epsilon;
	printf("Vector norm: %f Expected: %f\n", nm, 16.1245);


	return true;
}

int main(int argc, char * argv[]){
	int p = 6;
	bool norm_test = test_norm(p);
	printf("Norm test passed: %d\n", norm_test);

	bool matrix_test = test_matrix(p);
	printf("Matrix test passed: %d\n", matrix_test);

	bool restrict_test = test_restrict(p);
	printf("Restriction test passed: %d\n", restrict_test);

	bool interpolate_test = test_interpolate(p);
	printf("Interpolation test passed: %d\n", interpolate_test);


	return 0;
}
